Discrete Marks

Dependencies

Code
source("utils.R")

Setup

Code
spe <- readRDS("../data/spe.rds")

Multitype point process

Code
sub <- spe[, spe$sample_id == "0.01"]
(pp <- .ppp(sub, marks = "cluster_id"))
Marked planar point pattern: 6111 points
marks are of storage type  'character'
window: rectangle = [1222.5635, 3012.4248] x [-3993.535, -2202.755] units

Multitype and Multivariate viewpoint

A pattern with multiple type of points, e.g. cell types, can be seen in different ways. One the one hand, the multitype approach assumes that the points \(x\) were recorded together with with their labels \(m\) and that they were generated at the same time. The locations and labels therefore have a joint distribution \(P(X,M)\). On the other hand one can assume that the pattern with multiple types of points is a combination of several distinct point patterns, one for each type of point. This is the multivariate approach and the different point patterns \(A\) and \(B\) form a joint distribution \(P(A,B)\). To test if the labels depend on the location one can assume the following null hypotheses (Baddeley, Rubak, and Turner 2015, 565–67):

  • complete spatial randomness and independence (CSRI): the points are distributed at random; the type of each points is randomly allocated; independence between points of different types; allocation of the types independently of the other points and of its location.
  • random labeling: each point is assigned a type at random independently of its location
  • independence of components: the points of different types are independent of each other.

Apart from CSRI is is also important for the analysis if we can assume stationarity, i.e. the distribution of the points is the same everywhere in the window.

For simplicity, we will focus on three cell types of our point pattern: Ependymal, OD Mature and Microglia. This is appropriate if we assume that the point processes are independent. We could also assume that they come from the same process. In this case we have to check the stationarity assumption of the pattern.

Code
marks(pp) <- factor(marks(pp))
selection <- c('OD Mature', 'Ependymal', 'Microglia')

pp_sel <-  subset(pp, marks %in% selection, drop = TRUE)
Code
pp_sel |> as.data.frame() |> 
  ggplot(aes(x = x, y = y, color = marks)) +
  geom_point() +
  theme_minimal() +
  coord_fixed() +
  scale_color_brewer(palette = "Set1")

The summary of pp (point pattern) object returns general properties, plus intensities, combined and per mark type.

Code
summary(pp)
Marked planar point pattern:  6111 points
Average intensity 0.001906561 points per square unit

Coordinates are given to 4 decimal places

Multitype:
            frequency  proportion    intensity
Ambiguous         773 0.126493200 2.411670e-04
Astrocyte         646 0.105711000 2.015445e-04
Endothelial       469 0.076746850 1.463225e-04
Ependymal         268 0.043855340 8.361288e-05
Excitatory       1180 0.193094400 3.681463e-04
Inhibitory       1965 0.321551300 6.130571e-04
Microglia          97 0.015873020 3.026287e-05
OD Immature       200 0.032727870 6.239767e-05
OD Mature         457 0.074783180 1.425787e-04
Pericytes          56 0.009163803 1.747135e-05

Window: rectangle = [1222.5635, 3012.4248] x [-3993.535, -2202.755] units
                    (1790 x 1791 units)
Window area = 3205250 square units

To get the overall intensity the individual intensities can be summed up. Assuming that the the multitype process is first order stationary (i.e. each sub-process is stationary) the individual intensities sum up to the intensity of the unmarked point process (Baddeley, Rubak, and Turner 2015, 574ff.).

Code
sum(intensity(pp)) == intensity(unmark(pp))
[1] TRUE

The stationarity assumption is not appropriate in all cases. To assess stationarity visually, we can plot the kernel density estimates per type.

Code
ppls <- split(pp_sel) # split by mark
plot(density(ppls))

Ependymal and OD Mature cells are cleary inhomogeneous, while for Microglia cells it is not so clear and we could assume homogeneity, especially as there are no cells in the bottom middle of the window.

To further inverstiagte the spatial arrangement of the different cell types we can calculate the relative risk, i.e., the probability of observing a given celltype at a given location. It is calculated using the function relrisk. The bandwidth for smoothing is calculated with bw.relrisk and might need to be adjusted (Baddeley, Rubak, and Turner 2015, 577–83).

The relrisk function further gives us the dominant mark for different regions of the tissue of interest. This could be interesting in the annotation of spatial domains. It indicates at each location, which cell type is most likely to occur.

Code
rpd <- relrisk(pp_sel, diggle = TRUE)
dom <- im.apply(rpd, which.max)
dom <- eval.im(factor(dom, levels = seq_along(levels(unique(marks(pp_sel)))),
                      labels = levels(unique(marks(pp_sel)))))
plot(dom,las=2,main="Dominant mark")

Correlation and spacing

Nearest neighbourhood contingency

To further investigate the spatial distribution of the marks we can investigate the nearest neighbourhood of each cell type. One possibility is to work with nearest neighborhood contingency tables developed by (Dixon 2002). The statistical tests are implemented in the R package dixon (Cruz 2008).

The measure of segregation \(S\) is defined in (Dixon 2002) as

\[S_{i,j}= \frac{\log[(N_{i,j}/(N_i−N_{i,j})]}{[(N_i−1)/(N−N_i)]}\] where \(N_i\) is the number of individuals \(i\), \(N_{i,j}\) is the number of individuals of type \(i\) with a nearest neighbor of type \(j\), and \(N\) is the total number of individuals.

A value of \(S=0\) is consistent with random labeling. A value larger than 0 indicates that the two types are more segregated than expected by chance, the larger the value the more segregated. Note that segregated means that it is more likely to expect a neigbour of type \(j\) than by chance. In the case that the neigbour is of the same type this is equivalent to “clustering”. On the other hand if \(S<0\) it indicates that type \(j\) is less likely to be a neigbour than by chance. The P-values are calculated using expected numbers of nearest neighbors under the null hypothesis of random labeling using a Monte-Carlo simulation and assumes an asymptotic \(\chi^2\) distribution.

Code
out <- dixon(as.data.frame(pp_sel), nsim = 99)
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.
Code
out$tablaZ %>% 
  arrange(desc(abs(`Z `))) %>%
  select(-`  p-val.Nobs`)
       From        To     Obs.Count     Exp. Count    S      Z    p-val.Z
1 Ependymal Ependymal           262          87.16  1.96  20.04    0.0000
2 Ependymal OD Mature             3         149.18 -2.04 -16.87    0.0000
3 OD Mature Ependymal             8         149.18 -1.43 -14.26    0.0000
4 OD Mature OD Mature           380         253.83  0.60  11.43    0.0000
5 Ependymal Microglia             3          31.66 -1.07  -5.66    0.0000
6 Microglia Ependymal             9          31.66 -0.68  -4.92    0.0000
7 Microglia OD Mature            67          53.99  0.25   2.60    0.0094
8 Microglia Microglia            21          11.34  0.32   2.50    0.0124
9 OD Mature Microglia            69          53.99  0.12   2.35    0.0190

In this table we see that most Ependymal cells are very clustered, while Microglia are more evenly distributed. Further we see that it is less likely to find a Ependymal cells next to a OD mature cells than by chance.

OD Mature cells show this interesting characteristic that they are clustered in some parts of the tissue and more evenly distributed in other parts of the structure. This characteristic is not visible in the table. The statistic also considers only the nearest neighbour and ignores neighbours that are further away.

Summary functions for pairs of types

Similar to the simple case without marks, it is possible to estimate summary functions. In particular, summary functions between different marks can be calculated. Note that the canonical functions assume that the multi-type process is stationary.

Cross K-function

The cross K-function is a summary function that measures the average number of points of type j within a distance r of a point of type i. The formula is given by:

\[ K(r) = \frac{1}{\lambda_j} \mathbb{E} [t(u,r,X^{j})|u \in X^{i}], \]

where \(X^{i}\) is the point pattern of type \(i\) and \(t(u,r,X^{j})\) is the number of points of type \(j\) in a circle of radius \(r\) around \(u\) (Baddeley, Rubak, and Turner 2015, 594–95).

First, we plot an overview over the cross K function for the different types.

Code
plotCrossAll <- function(ppp, fun, edgecorr){
  nMarks <- length(unique(marks(ppp)))
  Fall <- alltypes(ppp, fun)
  
  # Create a list of ggplot objects using lapply
  plot_list <- lapply(Fall[["fns"]], function(res) {
    ggplot(res, aes(x = r, y = .data[[edgecorr]])) +
      geom_line(linewidth = 1) +
      geom_line(aes(x = r, y = theo), 
                linetype = "dotted", linewidth = 1) +
      geom_line() +
      labs(title = attributes(res)$yexp) +
      theme_minimal()
  })
  
  p <- wrap_plots(plot_list, ncol = nMarks) + 
    plot_layout(guides = "collect") & theme(legend.position='bottom')
  return(p)
}
Code
plotCrossAll(pp_sel, "Kcross.inhom", "iso")

The diagonal of the cross K-function plot shows the K-function for the different marks (indication of Poisson or non-Poisson point processes). Off-diagonal panels give indication of independence of points when the number of points follows the expected K-function but does not imply that the individual marks follow a Poisson process. If the types are independent, they are also uncorrelated.

It is important to remember that the cross K-function assumes that the multitype process is stationary. If this is not the case, there is a risk in misinterpreting the results. The problem is the confounding between clustering and inhomogeneity, c.f. (Baddeley, Rubak, and Turner 2015, 151–52) . In the example above, assuming that the process is inhomogeneous, the Ependymals cells appear to be regularly spaced, which seems counter intuitive. However, this is the result of the pattern being inhomogneous with spatially varying intensity. When accounting for this, the pattern is more regular than expected under an inhomogeneous point process. The estimation of the inhomogeneous cross functions is not straightforward and results change based on the estimation of the local intensity and the edge correction, c.f. (Baddeley, Rubak, and Turner 2015, 605).

In this overview, we can see that there is indication that Microglia and OD Mature cells are independent of each other. The other types seem to be dependent on each other. Let’s focus a bit more on the relationship between Ependymal and the other two cell types. We will also calculate confidence intervals for the different cross K-functions. We have already seen that our dataset most likely does not satisfy the assumption of stationarity. For this reason, we will calculate further calculate the inhomogeneous cross K-function.

Code
plotCrossMetric <- function(ppp, fun, from, to, edgecorr){
  lce <- lohboot(ppp, fun, from = from, to = to)
  p <- ggplot(lce, aes(x = r, y = .data[[edgecorr]])) +
    geom_line(size = 1) +
    geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.25)+
    geom_line(aes(x = r, y = theo), linetype = "dotted", size = 1) +
    geom_line() +
    labs(title = attributes(lce)$yexp) +
    theme_minimal()
  return(p)
}

p_epen_od <- plotCrossMetric(pp_sel, "Kcross.inhom", 
                             "Ependymal", "OD Mature", "iso")
p_epend_micro <- plotCrossMetric(pp_sel, "Kcross.inhom", 
                                 "Ependymal", "Microglia", "iso")
Code
# fig-width: 10
# fig-height: 10
p_epen_od + p_epend_micro

Remember that the dashed line represents the assumption of a multitype Poisson process. If the line lies above the dotted line, there is indication of clustering while if the line is below the dotted line there is indication of repulsion. In the plot above we can see that there is indication of clustering between Ependymal and OD Mature cells while there is indication of repulsion between Ependymal and Microglia cells.

Cross L-function

Alternatively the L cross function with similar interpretation can be calculated using the Lcross function (Baddeley, Rubak, and Turner 2015, 596ff).

Code
# fig-width: 10
# fig-height: 10
p_epen_od + p_epend_micro

Mark connection function

The mark connection function is the cross pair-correlation function, i.e. the generalization of the pair correlation function to a multitype point processes, divided by the unmarked pair-correlation function. It can be interpreted as the conditional probability that two points a distance \(r\)apart have labels of type A and of type B, given the presence of those points (Baddeley, Rubak, and Turner 2015, 596–97).

Code
plotCrossAll(pp_sel, "markconnect", "iso") + 
  scale_y_continuous(limits = c(0, 1))

The dashed lines indicate expected values under random labeling. The values measure dependence (or association) between the different labelled points. Positive values indicate that nearby points are more likely to have different types than expected by chance. This positive association between different cell types does not necessarily imply dependence, as it could be influenced by a negative association between cells of the same type, as it it could be the case for the Microglia cells. Furthermore, as the calculation is based on the \(K\) function, the mark connection function assumes homogenity.

Cross F-function (empty space function), cross G-function (Nearest-neighbor function) and cross J-function

The cross F-function is the cumulative distribution function of the distance from a location to the nearest point of the same type. For each type \(i\), it is defined as:

\[F_i(r) = \mathbb{P}\{d(u,X^{i}\leq r\}.\]

The cross G-function is the cumulative distribution function of the distance from a location to the nearest point of another type and is defined as:

\[G_{ij}(r) = \mathbb{P}\{d(x,X^{(j)} \setminus u \leq r \mid X^{(i)} \ \text{has a point at u}).\]

If the points are independent of each other, the G and F function are identical. Both assume that the process is stationary.

There exists a difference in the interpretation of the theoretical values of the K-cross and the G-cross function. For the K-cross, the theoretical value indicates independence between marks while for the G-cross the theoretical value is consistent with the assumption that the points of type j are Poisson in addition to being independent of the points of type \(i\) (Baddeley, Rubak, and Turner 2015, 597 ff).

The cross J-function is defined as:

\[J_{ij}(r) = \frac{1-G_{ij}(r)}{1-F_{j}(r)}\]

and summarizes the interpoint dependence between type \(i\) and \(j\). Under the hypothesis of independent components, i.e., that the point processes of each type are independent, the G-function is equivalent to the F-function and the J-function is equal to 1 (Baddeley, Rubak, and Turner 2015, 597 ff).

Dot functions

For each K-, G- and J- function, there also exist dot functions, which measure distances from points of one type to points of any type. These functions allow us to measure the dependence of one mark with all other marks at once. For expample, the K-dot function represents the expected number of an other point within distance \(r\) of a typical point of type \(i\) (Baddeley, Rubak, and Turner 2015, 600 ff).

Code
plotCrossAll(pp_sel, "Kdot.inhom", "iso")

The dot functions are useful summary statistics to analyse the dependence of one mark with all other marks.

Summary function within and between types

In our original dataset, we have a large number of different marks. We picked three OD mature, Ependymal and Microglia for illustrative purposes. An alternative to looking at all cross summary function combinations, it is possible to compare between and within types (Baddeley, Rubak, and Turner 2015).

Mark equality function

The Mark or Type Equality function for a stationary multitype point process measures the correlation between types of two points separated by distance r. It is the sum of the mark connection function of all pairs of points of the same type.

If k < 1, points at distance r are less likely than expected to be of the same type. If > 1, they are more likely to be of the same type. The value 1 indicates a lack of correlation (Baddeley, Rubak, and Turner 2015, 603 ff).

Code
plotMarkCorr <- function(pp, edgecorr = "iso") {
    me <- markcorr(pp)
    ggplot(me, aes(x = r, y = .data[[edgecorr]])) +
        geom_line(size = 1) +
        geom_line(aes(x = r, y = theo), linetype = "dotted", size = 1) +
        geom_line() +
        labs(title = attributes(me)$yexp) +
        theme_minimal()
}

plotMarkCorr(pp_sel)

We can see that in our dataset that it the more likely it is to find points of the same type at shorter distances. The curve never crosses the dashed line at 1, which means that it is generally more likely to find points of the same type at any distance than expected by chance.

Tests of randomness and independence

In a multitype point process, there are usually two interesting hypotheses:

  • random-labeling hypothesis: the allocation of labels to the points is random
  • independent component hypothesis: there is independence between different type of points

If both statments are correct, the point pattern is considered to be complete spatially random and independent (CSRI), the marked analog to complete spatial randomness (CSR) (Baddeley, Rubak, and Turner 2015, 605 ff).

Testing random labelling

The random labeling test is most logical when the marks represents its status, which is not most appropriate assumption when considering cell types. Testing for random labeling can be done using permutation test, in which the labels are randomly permuted. Random labeling can be assumed if the permuted datasets are statistically equivalent to the original dataset (Baddeley, Rubak, and Turner 2015, 609 ff).

Testing the indepenence of components assumption

The i-to-j functions are useful to test the independence of different subprocesses. If the processes of type i and j are independent, then \(K_{ij} = \pi r^2, G_{ij}(r) = F_{j}(r), J_{ij}(r) \equiv 1\). Alternatively, randomization tests can be used in which simulated patterns from the dataset are generated and randomly split into subpatterns. These are then compared to the null hypothesis in which all subpatterns should be statistically equivalent to the original. However, this approach assumes stationarity and there is a need to handle edge effects (Baddeley, Rubak, and Turner 2015, 606 ff).

Code
plotEnvCross <- function(pp, i, j, fun, nsim = 39, radius = 150, global = FALSE){
  pp_scaled <- rescale(pp)
  E1 <- envelope(pp_scaled, fun, nsim=nsim, i=i, j=j,
                 simulate=expression(rshift(pp_scaled, radius = radius)), global = global)
  p <- ggplot(E1, aes(x = r, y = .data[["mmean"]])) +
    geom_line(size = 1) +
    geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.25)+
    geom_line(aes(x = r, y = obs), linetype = "dotted", size = 1) +
    geom_line() +
    labs(title = attributes(E1)$yexp) +
    theme_minimal()
  return(p)
}

pEnv <- plotEnvCross(pp_sel, fun = "Kcross.inhom", 
                     "Ependymal", "OD Mature", nsim = 39, radius = 150)
Code
pEnv

Code
plotEnvCross(pp_sel, fun = "Kcross.inhom", 
             "Ependymal", "OD Mature", nsim = 39, radius = 150, global = TRUE)
Generating 78 simulations by evaluating expression (39 to estimate the mean and 
39 to calculate envelopes) ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 
78.

Done.

We have indication that the independence-of-components assumption should not be rejected. Therefore, we conclude that Ependymal and OD Mature cells are independent.

Assuming stationarity of the total pattern

As outlined above, the cross correlation and spacing functions assume stationarity, which is not given in our dataset, when we look at the different patterns individually. However, using the total (unmarked) pattern, we could assume stationarity.

Code
plot(density(unmark(pp)))

Let’s look at the homogeneous cross K-function.

Code
plotCrossAll(pp_sel, "Kcross", "iso")

The result is different from the previous analysis. The Ependymal cells now appear to be clustered. This is because stationarity assumes that Ependymal cells could theoretically be present in the total observation window. If this assumption is justified, depends on the context and research question. The interpretation of the results should always be done with the assumption of stationarity or inhomogeneity in mind and should be reported in an analysis.

Session info

Code
sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Zurich
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] magrittr_2.0.3                 stringr_1.5.0                 
 [3] dixon_0.0-8                    splancs_2.01-44               
 [5] spdep_1.2-8                    spData_2.3.0                  
 [7] tmap_3.3-4                     scater_1.28.0                 
 [9] scran_1.28.2                   scuttle_1.10.3                
[11] SFEData_1.2.0                  SpatialFeatureExperiment_1.2.3
[13] Voyager_1.2.7                  rgeoda_0.0.10-4               
[15] digest_0.6.33                  ncf_1.3-2                     
[17] sf_1.0-14                      reshape2_1.4.4                
[19] patchwork_1.1.3                STexampleData_1.8.0           
[21] ExperimentHub_2.8.1            AnnotationHub_3.8.0           
[23] BiocFileCache_2.8.0            dbplyr_2.3.4                  
[25] RANN_2.6.1                     seg_0.5-7                     
[27] sp_2.1-1                       rlang_1.1.1                   
[29] ggplot2_3.4.4                  dplyr_1.1.3                   
[31] mixR_0.2.0                     spatstat_3.0-6                
[33] spatstat.linnet_3.1-1          spatstat.model_3.2-6          
[35] rpart_4.1.19                   spatstat.explore_3.2-3        
[37] nlme_3.1-162                   spatstat.random_3.1-6         
[39] spatstat.geom_3.2-5            spatstat.data_3.0-1           
[41] SpatialExperiment_1.10.0       SingleCellExperiment_1.22.0   
[43] SummarizedExperiment_1.30.2    Biobase_2.60.0                
[45] GenomicRanges_1.52.1           GenomeInfoDb_1.36.4           
[47] IRanges_2.34.1                 S4Vectors_0.38.2              
[49] BiocGenerics_0.46.0            MatrixGenerics_1.12.3         
[51] matrixStats_1.0.0             

loaded via a namespace (and not attached):
  [1] spatstat.sparse_3.0-2         bitops_1.0-7                 
  [3] httr_1.4.7                    RColorBrewer_1.1-3           
  [5] tools_4.3.1                   utf8_1.2.3                   
  [7] R6_2.5.1                      HDF5Array_1.28.1             
  [9] mgcv_1.8-42                   rhdf5filters_1.12.1          
 [11] withr_2.5.1                   gridExtra_2.3                
 [13] leaflet_2.2.0                 leafem_0.2.3                 
 [15] cli_3.6.1                     labeling_0.4.3               
 [17] proxy_0.4-27                  R.utils_2.12.2               
 [19] dichromat_2.0-0.1             scico_1.5.0                  
 [21] limma_3.56.2                  rstudioapi_0.15.0            
 [23] RSQLite_2.3.1                 generics_0.1.3               
 [25] crosstalk_1.2.0               Matrix_1.5-4.1               
 [27] ggbeeswarm_0.7.2              fansi_1.0.5                  
 [29] abind_1.4-5                   R.methodsS3_1.8.2            
 [31] terra_1.7-55                  lifecycle_1.0.3              
 [33] yaml_2.3.7                    edgeR_3.42.4                 
 [35] rhdf5_2.44.0                  tmaptools_3.1-1              
 [37] grid_4.3.1                    blob_1.2.4                   
 [39] promises_1.2.1                dqrng_0.3.1                  
 [41] crayon_1.5.2                  lattice_0.21-8               
 [43] beachmat_2.16.0               KEGGREST_1.40.1              
 [45] magick_2.8.0                  pillar_1.9.0                 
 [47] knitr_1.44                    metapod_1.7.0                
 [49] rjson_0.2.21                  boot_1.3-28.1                
 [51] codetools_0.2-19              wk_0.8.0                     
 [53] glue_1.6.2                    vctrs_0.6.4                  
 [55] png_0.1-8                     gtable_0.3.4                 
 [57] cachem_1.0.8                  xfun_0.40                    
 [59] S4Arrays_1.0.6                mime_0.12                    
 [61] DropletUtils_1.20.0           units_0.8-4                  
 [63] statmod_1.5.0                 bluster_1.10.0               
 [65] interactiveDisplayBase_1.38.0 ellipsis_0.3.2               
 [67] bit64_4.0.5                   filelock_1.0.2               
 [69] irlba_2.3.5.1                 vipor_0.4.5                  
 [71] KernSmooth_2.23-21            colorspace_2.1-0             
 [73] DBI_1.1.3                     raster_3.6-26                
 [75] tidyselect_1.2.0              bit_4.0.5                    
 [77] compiler_4.3.1                curl_5.1.0                   
 [79] BiocNeighbors_1.18.0          DelayedArray_0.26.7          
 [81] scales_1.3.0                  classInt_0.4-10              
 [83] rappdirs_0.3.3                goftest_1.2-3                
 [85] fftwtools_0.9-11              spatstat.utils_3.0-3         
 [87] rmarkdown_2.25                XVector_0.40.0               
 [89] htmltools_0.5.6.1             pkgconfig_2.0.3              
 [91] base64enc_0.1-3               sparseMatrixStats_1.12.2     
 [93] fastmap_1.1.1                 htmlwidgets_1.6.2            
 [95] shiny_1.7.5.1                 DelayedMatrixStats_1.22.6    
 [97] farver_2.1.1                  jsonlite_1.8.7               
 [99] BiocParallel_1.34.2           R.oo_1.25.0                  
[101] BiocSingular_1.16.0           RCurl_1.98-1.12              
[103] GenomeInfoDbData_1.2.10       s2_1.1.4                     
[105] Rhdf5lib_1.22.1               munsell_0.5.0                
[107] Rcpp_1.0.11                   ggnewscale_0.4.9             
[109] viridis_0.6.4                 stringi_1.7.12               
[111] leafsync_0.1.0                zlibbioc_1.46.0              
[113] plyr_1.8.9                    parallel_4.3.1               
[115] ggrepel_0.9.4                 deldir_1.0-9                 
[117] Biostrings_2.68.1             stars_0.6-4                  
[119] splines_4.3.1                 tensor_1.5                   
[121] locfit_1.5-9.8                igraph_1.5.1                 
[123] ScaledMatrix_1.8.1            BiocVersion_3.17.1           
[125] XML_3.99-0.14                 evaluate_0.22                
[127] BiocManager_1.30.22           httpuv_1.6.11                
[129] purrr_1.0.2                   polyclip_1.10-6              
[131] rsvd_1.0.5                    lwgeom_0.2-13                
[133] xtable_1.8-4                  e1071_1.7-13                 
[135] RSpectra_0.16-1               later_1.3.1                  
[137] viridisLite_0.4.2             class_7.3-22                 
[139] tibble_3.2.1                  memoise_2.0.1                
[141] beeswarm_0.4.0                AnnotationDbi_1.62.2         
[143] cluster_2.1.4                

References

Baddeley, Adrian, Ege Rubak, and Rolf Turner. 2015. Spatial Point Patterns. 1st ed. CRC Interdisciplinary Statistics Series. CRC Press, Taylor & Francis Group.
Cruz, Marcelino de la. 2008. “Métodos Para Analizar Datos Puntuales.” In Introduccion Al Analisis Espacial de Datos En Ecologia y Ciencias Ambientales: Metodos y Aplicaciones, 76–127.
Dixon, Philip M. 2002. “Nearest-Neighbor Contingency Table Analysis of Spatial Segregation for Several Species.” Écoscience 9 (2): 142–51. https://doi.org/10.1080/11956860.2002.11682700.